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The Bethe-Salpeter equation is a widely used approach to describe optical excitations in bulk 
semiconductors. It leads to spectra that are in very good agreement with experiment, but the price 
to pay for such accuracy is a very high computational burden. One of the main bottlenecks is the 
large number of fc-points required to obtain converged spectra. In order to circumvent this problem 
we propose a strategy to solve the Bethe-Salpeter equation based on a double-grid technique coupled 
to a Wannier interpolation of the Kohn-Sham band structure. This strategy is then benchmarked 
for a particularly difficult case, the calculation of the absorption spectrum of GaAs, and for the well 
studied case of Si. The considerable gains observed in these cases fully validate our approach, and 
open the way for the application of the Bethe-Salpeter equation to large and complex systems. 



PACS numbers: 78.20.-c Optical properties of bulk materials and thin films; 78.20.Bh Theory, models, and 
numerical simulation: 



Optical spectra are important for the characterization 
and prediction of material properties, as optical exci- 
tations are at the core of e.g., light-emitting devices, 
laser technology, and photovoltaics. For extended sys- 
tems many-body perturbation theory^ a Green's func- 
tion based approach, is the most accurate method to cal- 
culate optical properties. Perhaps inevitably, it is also 
one of the most computationally costly methods avail- 
able to the community. It involves the solution of an 
equation of motion for the two-particle Green's function, 
the Bethe-Salpeter equation (BSE), that describes cou- 
pled and correlated electron-hole excitations^ 

The standard numerical techniques used to solve the 
BSE are based on an expansion of the relevant quanti- 
ties in electron-hole states (needing therefore both filled 
and empty states), and require a very dense fc-point sam- 
pling of the Brillouin zone (BZ). Typically, the number 
of electron-hole states used in the expansion can be rela- 
tively small if one is only interested in the visible spectra, 
but the number of /c-points can easily reach several thou- 
sands. Some approaches have been put forward to reduce 
the computational burden of the BSE. For example, the 
number of k-points can be reduced by interpolating the 
interaction integrals in fc-space^ while recent implemen- 
tations allow for the complete exclusion of empty states^ 

It is well known that optical spectra are very sensitive 
to the fc-point sampling^— A common approach to alle- 
viate the problem is the use of arbitrarily shifted fc-point 
grids, that often yield sufficient sampling of the Brillouin 
zone while keeping the number of fc-points manageable. 
Such a shifted grid, indeed, does not use the symmetries 
of the Brillouin zone and guarantees a maximum number 
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of non-equivalent fc-points accelerating spectrum conver- 
gence^ However, it might induce artificial splitting of 
normally degenerate states, and thus produce artifacts 
in the spectrum^ such as the splitting of some peaks 
or even the appearance of spurious excitations in some 
directions. Of course, these artifacts (slowly) disappear 
with increasing density of fc-points^^ and consequent in- 
crease of the computational burden. In view of that a 
very dense fe-point sampling is crucial to obtain an accu- 
rate lineshape, including the correct peak positions^^ 
but very hard to achieve in practical calculations. 

In this Article, we propose a new strategy to solve the 
BSE equation that alleviates the need for dense fc-point 
grids. The independent-particle part of the BSE is first 
evaluated on a very dense fc-grid (40x40x40 in the exam- 
ple below) by making use of Wannier interpolation.— The 
BSE is then solved in a unshifted coarse fc-grid (10x10x10 
in the example below) using a double-grid technique to 
take into account the fast changing independent-particle 
contribution. This approach is simple to implement, and 
leads to a considerable gain in computational time. 

In the following, we start by presenting a short review 
of the theoretical ingredients for the description of optical 
spectra within the BSE^i2 We then discuss our approach 
and prove its usefulness with a notoriously difficult exam- 
ple, the calculation of the optical absorption spectrum of 
the standard semiconductor GaAs. The subsequent dis- 
cussion of bulk silicon concludes the benchmark. 

The optical absorption spectrum is described by the 
imaginary part of the macroscopic dielectric function 
£m(w) in the long wavelength limit, which in turn can be 
obtained from the two-point contraction of the reducible 
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four-point polarizability L, 



In this basis the polarizability is written as 



e M M = 1-lim v(q)X [ drdr' e- lq{r - r,) L(r,r,r' ,r';u), 
<?^o J 

(1) 

with the Coulomb potential v = A-k/(G + q) 2 , the trans- 
ferred momentum q, and A the direction of light polar- 
ization. The quantity L satisfies the BSE, a Dyson like 
equation, 

L(1,2,3,4) = L°(1,2,3,4)+ 

<f(5678)L°(l, 2, 5, 6) S(5, 6, 7, 8) L(7, 8, 3, 4), (2) 



with the abbreviation of space, spin and time coor- 
dinates (1) = (t*i, <7i, ti), and where L°(l,2, 3,4) = 
iG(l, 3)G(4, 2) is the independent particle polarizability, 
expressed as a product of single-particle Green's func- 
tions. Equation @ describes the effects of the electron- 
hole interaction mediated by the BSE kernel 3 = v — W 
that is composed of (i) a bare, repulsive short-range ex- 
change term, that includes the microscopic components 
of the Coulomb interaction, i.e. = vg',vg=q = 0, 

and (ii) an attractive, static screened Coulomb potential 
W, the direct term, arising from the variation of the self- 
energy. Dynamical effects due to the self-energy influence 
both the (single) quasiparticle renormalization and the 
excitonic two-body interaction Wi 13 ' 14 In response cal- 
culations for semiconductors they partially cancel each 
other, which justifies the commonly employed approxi- 
mation of a static W and neglected quasiparticle renor- 
malization, but in general this is not true for metalsi^ 4 - 

As the interaction is instantaneous, only two time- 
variables of the initial four remain. And, due to the 
time-translation invariance, L and L° depend only on 
the relative time-difference. A time-energy Fourier trans- 
formation then turns the polarizability into a function 
of a single frequency L( 1,2,3, 4; a;) where from now on 
(1) = (ri,<7i) only. 

By taking advantage of the two-particle nature of the 
BSE, all expressions are conveniently written on the ba- 
sis of the electron- hole vertical transition space composed 
of N v valence bands, N c conduction bands, and N^ z k- 
points in the whole BZ. The dimension of this basis is 
4xN v x N c x N^ z . In the case of vanishing spin-orbit cou- 
pling, the BSE can be separated into the two subspaces of 
singlet and triplet excitons^ each of them with dimension 
2xN v xN c xN^ z . The basis sets that span these subspaces 
are constructed from pairs of single-particle states 4> n ^k 
with n as band index and k as fc-point and spin variable, 
such that 



$jc(ri, r 2 ) = 4> c ,k(ri) ■ <t>* k (r 2 ), 



(3) 



with the short hand notation K = (c,v,k), and c and 
v running over indices of conduction and valence bands, 
respectively 



Lk!,k 2 (u)= J dr 1 dr 2 dr 3 dr 4 

$* Kl (n, r 2 )L(l, 2, 3, 4; w)$x 2 (t- 3 , r 4 ), (4) 
and L°, that is now diagonal, reads 
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where / denotes the occupation number. The infinitesi- 
mal r) shifts the pole u> = e ci k t — £ vi ki away from the real 
axis, and is thus responsible for a finite life-time of the 
excitation. 

In this basis, the BSE becomes a matrix equation 



L%c u K 3 -K 3 ,K i L KifK2 



(6) 



where we used Einstein's notation for summations over 
repeated indices in the tensor products, and omitted the 
explicit energy dependence for clarity. In the following 
we adopt the notation of O for the matrix representation 
in the electron-hole basis of an arbitrary operator O. 

It can be shown that the kernel S couples pairs of ex- 
citations (vc) with (v'c 1 ), but also (vc) with (cV), lead- 
ing to the so-called resonant and coupling terms, respec- 
tively. 4 Here 1 we mala 1 use of the so-called Tamm-Dancoff 
approximation, and neglect the latter. We thus arrive at 
a Hilbert space of dimension N v xN c x N^ z , regardless of 
the symmetries of the /c-grid. Note that these are stan- 
dard approximations for the solution of the BSE. 

Equation ([6]) can be solved symbolically, yielding for 
each frequency u> 



J Ki,K 2 



L 
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K!,K 3 K3,K 2 



(7) 



To circumvent the inversion of 1 — L S , the usual proce- 
dure is to rewrite the matrix into a two-particle Hamilto- 
nian whose diagonalization gives the excitonic eigensys- 
tem used to express the polarizability for all frequencies 
at once. 

As it will become clear in the following, we stick to 
the inversion scheme by taking advantage of the series 
expansion of Eq. (J7J) , 



M, (8) 



that is interrupted at convergence. If, however, a conver- 
gence is not attained, i.e. the assumption of the expand- 
ability of Eq. ([7]) is falsified a posteriori, we perform the 
full inversion. 

The solution of Eq. © has two distinct bottlenecks in 
terms of computational cost. First, the calculation of ma- 
trix S is very time consuming. Second, its storage needs 
large quantities of memory. In view of that, reducing the 
number of Appoints is a major issue. To this end, Rohlf- 
ing and Louie^ employed a double-grid technique where 
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the kernel 3 is calculated on a coarse grid with its sub- 
sequent interpolation onto a fine fc-point mesh where the 
BSE is solved. This approach helps reducing the time 
necessary to compute 3, but it is less helpful to save 
memory, since it requires the storage of the computed 
and interpolated matrix elements of the kernel. Its use is 
justified with the authors' observation that 3 varies little 
with respect to the fc-points, as the single-particle wave- 
functions cj) n k are quite robust with respect to k (with the 
possible exception of sudden band crossings). In a sim- 
ilar spirit it has been prove d 15 ! 16 for the random phase 
approximation (RPA), that LP_ is a rapidly varying quan- 
tity and that could be correctly evaluated by performing 
additional Monte Carlo integrations on a large number of 
random k points. Note that LP_ can in principle be easily 
calculated for a large number of fc-points and bands. 

By taking into account these observations, we define 
in our approach two grids: (i) a coarse one with points k 
on which we calculate and store 3 and solve Eq. ([5]) , and 
(ii) a fine fc-grid with vectors k on which we compute 
L° . The mapping of L ^ ^ to L%c t k 2 is performed 

through a double-grid technique^ with a suitably chosen 
interpolation for the kernel. To simplify our approach we 
use the simplest zeroth order interpolation, that leads 
to averaging the finely resolved LP_ in a neighborhood 
around each point of the coarse grid. Consequently, this 
technique is expected to work if the oscillator strengths 
and 3 are smoothly varying functions of k. In practice, 
we define: 



L°(lo)k 1 ,k 2 
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where is the number of fc-points of the fine grid in the 
domain T>k around k of the coarse grid. We would like 
to note that we are averaging the polarization LP_ that 
has poles at the particle excitation energies, which is not 
equivalent to averaging the particle excitation energies 
that appear in the diagonal of the excitonic Hamiltonian. 
An arbitrary fc-point resolution of LP_ is possible once the 
respective single-particle energies e n ^ are available. 

In general the calculation of quasiparticle states on 
the k grid are not practical. Fortunately, there is a so- 
lution to this problem that relies on the interpolation 
of the (quasiparticle) electronic structure to a dense k- 
grid using maximally localized Wannier functions^ In 
this method, e £ at an arbitrary fc-point k is the result 
of (i) a rotation of the initial quasiparticle Hamiltonian 
into the Wannier basis, (ii) its Fourier interpolation to 
the fine grid of fc-points, and (iii) the diagonalization 
of the resulting Hamiltonian^ Note that, even if this 
procedure leads to an expression that has the form of a 
Slater-Koster tight-binding interpolation^ the obtained 
single-particle energies are calculated directly from the 
underlying ab initio eigensystem rather than just fitted 
to it. 

To illustrate the efficiency of our scheme we calculated 
the optical spectra of semiconducting GaAs, known for 
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FIG. 1. (Color online) Calculated RPA absorption spectra 
starting from GW corrected bands for bulk GaAs with N v = 
N c = 2 for various fc-grids, without (a) and with (b) double- 
grid method. The shaded area in (b) shows the spectrum 
corresponding to the 40x40x40 fc-grid of panel (a) for better 
comparison between the double-grid method and the standard 
calculations. 
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FIG. 2. (Color online) Calculated BSE absorption spectra 
for bulk GaAs with N v = 2 and N c = 3 using several k- 
grids without (a) and with (b) double-grid method. In (a) 
one spectrum is calculated on a symmetric 10x10x10 fc-grid, 
while the others are on shifted grids to accelerate convergence. 
The shaded area in both panels indicate the experimental 
spectrum at 22 K.— 



its slow convergence with respect to the fc-point sam- 
pling.— The Kohn-Sham band structure and wave- 
functions were obtained with density functional the- 
ory (DFT) within the local density approximation using 
norm conserving pseudopotentials with an energy cutoff 
of 14 Hartree and the experimental lattice constant of 
10.68 Bohr^i 
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For the DFT part we utilized the code ABINIT^ To 
obtain the energy bands [used in Eq. ©] we employ the 
code wannier9022 to perform a Wannier interpolation 
to a 40 x 40 x 40 regular fc-point grid in the whole BZ. 
Finally, optical spectra were calculated using the code 
Yambo— that uses the DFT Kohn-Sham wavefunctions 
and the interpolated single-particle energies as input. 

It is well known3^ that self-energy corrections in GaAs 
can be simulated by a rigid shift of the conduction bands. 
We therefore applied a scissor operator of 0.9 eV, that 
yields an overall agreement of the band dispersions within 
0.1 eV with the GW corrected bands, and a close agree- 
ment with experimental data as welli^— 

For all spectra in Fig. [T] and [2] we included the two 
highest valence bands and the two (three for BSE) lowest 
conduction bands, considering only the resonant part of 
the BSE kernel S, and used a Lorcntzian broadening of 
0.1 eV. Furthermore, we neglected spin-orbit coupling. 
Omitting local field effects (LFE), the non-interacting 
RPA spectra in Fig. [lja) arc obtained on symmetric 
Monkhorst-Pack (MP) grids£L With increasing fc-point 
resolution the spectrum converges to two main peaks at 
3.3 eV and 5.3 eV. By using our double-grid method, 
shown in panel (b), a 12 x 12 x 12 symmetric grid yields 
an equally well converged spectrum. We observe an ex- 
cellent agreement between the latter and the RPA done 
on a 40 x 40 x 40 grid (indicated by the shaded area) . 

In general, if one calculates independent-electron tran- 
sitions starting from GW corrected bands the oscilla- 
tor strength of the absorption spectrum is moved too 
high in energy compared to the experiment. The at- 
tractive net electron-hole interaction decreases the energy 
of the excited states and transfers oscillator strength to 
lower energies. This can be seen by comparing the non- 
interacting and interacting results of Figs. [1] and [3J re- 
spectively. 

Figure HJa) illustrates that a symmetric 10 x 10 x 10 
grid alone does not provide enough independent sam- 
pling points for a converged BSE spectrum. Shifting this 
grid in a direction different from the high symmetry di- 
rections provides 1 000 instead of only 47 nonequivalent 
sampling points (see Table H| . This leads to a spectrum 
that is sufficiently compatible with experimental results 
of Ref. [l^. Nevertheless, the low-energy region (peak 
at 1.9 eV) and the region between the two main transi- 
tions at 3.2 eV and 5.1 eV are still expected to change 
on a denser fc-grid. With our double-grid technique the 
BSE spectrum is converged even on a symmetric grid of 
10x10x10 [see Fig.[2](b)]. Contrary to the spectrum on an 
equally dense, but shifted grid, our spectrum is smooth 
in between the main transitions at 3.2 eV and 5.1 cV. It 
is noteworthy to mention that, for our scheme, a shifted, 
coarse fc-grid does not improve the convergence of the 
spectrum of GaAs. 

In the same fashion we calculated the RPA and BSE 
absorption spectra of Si, shown in Fig. EI To converge the 
ground state KS energies and wavefunctions with DFT 
we used an energy cutoff of 15 Ha and a lattice constant 
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FIG. 3. (Color online) Calculated RPA (a) and BSE (b) ab- 
sorption spectra with LFE starting from GW corrected bands 
for bulk Si with and without double-grid method employing a 
broadening of 0.05 eV. For (a) N v = N c = 2 were used, while 
for (b) N v = 2 and N c = 3. With the double-grid method 
both spectra (solid lines) converge faster. The shaded area in 
(b) shows the experimental spectrum.— 



of 10.2 Bohr obtained by crystal relaxation^ Similarly 
to GaAs, we took advantage of the good approximation 
of self-energy corrections by a rigid shift of the conduc- 
tion bands of 0.8 eV£ For the spectra we employed the 
two highest valence bands and the two (three for BSE) 
lowest conduction bands, and we included LFE by a di- 
electric matrix of size 51x51. Furthermore, we diminished 
the Lorentzian broadening to 0.05 eV in order to better 
resolve the first peak of the BSE spectrum at 3.5 eV. 

Consequently a high fc-point resolution of 60x60x60 was 
necessary to converge the RPA spectrum on a symmetric 
grid without the double-grid method. Figure [3][a) illus- 
trates the advantage of the double-grid technique as it 
warrants a converged RPA spectrum on a 12x12x12 grid. 
With this method also the BSE spectrum in Fig. 0(b) 
converged fast on a 10 x 10 x 10 grid, while the standard 
method of diagonalizing the BSE Hamiltonian is still far 
off convergence on the same fc-grid. Additionally, we ob- 
serve a good agreement of the converged BSE spectrum 
with experiment^ (indicated by the shaded area). 

A final remark for Si on the sampling of the dense fc- 
points k in Eq. @ is in order now. Instead of using a 
40x40x40 regular fc-grid of the full BZ (as for GaAs), 
we found it favorable to resort to a set of four 40 x 40 x 40 
shifted MP grids that respects the face-centered cubic 
symmetry of the Si crystal. Note that the calculations 
of the spectra using the double-grid technique were then 
again performed on unshifted, symmetric MP grids. 

Although a MP grid that respects the symmetries of 
the BZ does not reduce the dimension of the BSE kernel, 
its use is still advantageous for three reasons. Firstly, 
there is no artificial splitting of degenerate states^ Sec- 
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TABLE I. Number of fc-points in the (ir)reducible BZ N^ )BZ 
that are used for the calculation of the spectra in Figs. [U -H 

calc. Figs. fc-point grid N™ z A| z 

RPA Ha) 60x60x60 5 216 216 000 

RPA [Ha) 40x40x40 1 661 64 000 

RPA, dgrid CUb), |3ja) 12x12x12 72 1 728 

BSE (Ha) 10x10x10 shifted 1 000 1 000 

BSE, dgrid |{bjj|b) 10x10x10 47 1 000 

ondly, no artificial crystal anisotropy is introduced that 
has to be compensated by averaging the computed spec- 
tra over the three spatial directions of light polariza- 
tion£i2i Finally, in the calculation of the exchange term 
of S, symmetries of the BZ can be exploited, which trans- 
lates in a strong reduction of computational time. In Ta- 
ble U we have summarized the number of fc-points with 
and without considering symmetries of the BZ. 

In conclusion, we presented a double-grid method to 



solve the BSE on a coarse fc-point grid, where the av- 
erage of the strongly varying, but easily obtainable, 
independent-particle polarization is used. Converged 
spectra are reached for relatively small symmetric fc-point 
grids. This allows for a considerably faster calculation of 
the BSE kernel. The single-particle energy bands in a 
dense fc-point grid, the basic ingredient of our method, 
are not calculated directly, but are obtained through 
Wannier interpolation of the electronic band-structure. 
As examples, we discussed the convergence of the absorp- 
tion spectra of GaAs and Si with respect of the number 
of fc-points. The speed-up is considerable, and opens the 
way for the solution of the BSE equation in large, com- 
plex systems. 
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